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Two relations, the virial relation Madm = Mk and the first law in the form SAIadm = i^SJ, 
should be satisfied by a solution and a sequence of solutions describing binary compact objects in 
quasiequilibrium circular orbits. Here, Madm, Mk, J, and are the ADM mass, Komar mass, 
angular momentum, and orbital angular velocity, respectively. 5 denotes an Eulerian variation. 
These two conditions restrict the allowed formulations that we may adopt. First, we derive relations 
between Madm and Mk and between SMadm and Q,SJ for general asymptotically flat spacetimes. 
Then, to obtain solutions that satisfy the virial relation and sequences of solutions that satisfy the 
first law at least approximately, we propose a formulation for computation of quasiequilibrium binary 
neutron stars in general relativity. In contrast to previous approaches in which a part of the Einstein 
, equation is solved, in the new formulation, the full Einstein equation is solved with maximal slicing 

and in a transverse gauge for the conformal three-metric. Helical symmetry is imposed in the near 
zone, while in the distant zone, a waveless condition is assumed. We expect the solutions obtained 
in this formulation to be excellent quasiequilibria as well as initial data for numerical simulations of 
' binary neutron star mergers. 
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\l . A detailed theoretical understanding of evolution of close binary neutron stars is one of the most important goals in 
■ general relativity, since they are promising sources of gravitational waves for laser interferometric gravitational wave 
O , detectors such as LIGO, TAMA, GEO, and VIRGO [1,2]. For the inspiral phase in which the orbital separation a 
O^' is much larger than the radius i? of a neutron star, the orbital velocity is much smaller than the speed of light, and 
[ finite-size effects of neutron stars may be ignored. Thus, a post Newtonian study together with the point particle 
approximation is appropriate [3]. For a/R ^ 4, however, the post Newtonian and point particle approximations 
^ \ break down and numerical study is required to take into account the effect of tidal deformation of each star and full 
• ^ ■ effects of general relativity. The procedure to be adopted for such close orbits up to the merger is (i) to compute a 
' quasiequilibrium circular orbit at a distant orbit with a ^ AR and a ^ lOAf for which ratio of a radial approaching 
velocity to the orbital one will be small (less than 1%) [4], and then (ii) to perform a numerical relativity simulation 
adopting a distant quasiequilibrium with a ^ lOM as the initial condition [5,6]. In this paper, we focus on the 
formulation for computation of the quasiequilibrium in circular orbits. 

So far, the quasiequilibrium states of binary neutron stars have been widely computed in the so-called conformal 
fiatness approximation (or Isenberg-Wilson-Mathews formalism) [7-1 1], in which the conformal three- metric is as- 
sumed to be flat. The solution in this formulation satisfies the constraint equations of general relativity and, hence, it 
is fully general relativistic for the initial value problem. However, it is only an approximate quasiequilibrium solution 
for compact binaries, since conformally nonflat parts of the three- metric are not vanishing for the quasiequilibrium 
binaries. Thus, such approximation produces a systematic error of magnitude ~ (M/a)"^ for the solution of quasiequi- 
librium configurations [4] . The systematic error is also included in gravitational waves computed from the conformal 
flat data of the quasiequilibrium binary [4,12] and in the numerical results of fully general relativistic simulations 
started from initial conditions of the conformally flat quasiequilibria [5]. 

Formulations for computation of binary compact objects in quasiequilibrium circular orbits with a conformally 
nonflat three-metric have been proposed by several authors (e.g., [13-18,4]). A promising approach to this problem is 
to assume a helical Killing symmetry for the spacetime [17,18]. In this case, however, the solution contains standing 
gravitational waves in the whole spacetime and the averaged energy density of gravitational waves falls off as 
where r denotes a radial coordinate, resulting in an asymptotically nonflat spacetime. Thus, the solution obtained in 
such formulation is not physical in the distant wave zone, although the solution in the near zone and in a local wave 
zone would describe a realistic spacetime of binary compact objects. 
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In this paper, wc consider general relativistic formulations for computation of the quasiequilibrium circular orbits 
assuming that the spacetime is asymptotically flat. First, we require that the following two conditions should be 
satisfied for a solution and a sequence of the solutions of quasiequilibrium states: 

• A quasiequilibrium solution that is stationary in the corotating frame should satisfy a virial relation associated 
with the equality 

Madm = Mk, (1.1) 

of the ADM and Komar masses defined in Sec. III. 

• Binary compact objects inspiral adiabatically as a result of gravitational wave emission, conserving baryon mass, 
entropy, and vorticity. Thus, along a sequence of quasiequilibrium solutions, the first law should be satisfied. 
Here, the first law is written in the form 

^Madm = ^SJ, (1.2) 

where SMabm and 5 J are infinitesimal differences of the ADM mass and angular momentum along a quasiequi- 
librium sequence, and O is an orbital angular velocity. 

These two conditions are likely to be satisfied for binary neutron stars in nature. Thus, we should adopt a formulation 
which provides a solution and a sequence of the solutions that satisfy two conditions at least approximately. 

Based on this motivation, in this paper, we first derive relations for the differences, Madm — Mk and 6Madm — i^SJ, 
in arbitrary asymptotic flat spacetimes. The condition that the differences vanish can be used to restrict formulations 
that we can adopt. Using these conditions, several possible candidates for the formulations emerge. Among them, we 
propose a formulation in which helical symmetry is imposed only in the near zone instead of in the whole spacetime. 
Specifically, wc impose a mixed condition; a helical symmetry condition in the near zone and a wavclcss condition 
in the distant zone. To fix the gauge, we adopt the maximal slicing condition and a transverse gauge condition for 
the upper component of the conformal three-metric. In this case, all components of the Einstein equation reduce to 
elliptic equations as in the post Newtonian approximation. This implies that no standing waves appear in the wave 
zone, although in the near zone, gravitational-wave-like components are present. We expect the solutions obtained in 
this formulation to be excellent quasiequilibria as well as initial data for numerical simulations of binary neutron star 
mergers. 

The paper is organized as follows. In Sec. II, we describe the basic equations for quasiequilibria. In Sec. Ill, we 
derive a relation between Madm and Mk for arbitrary formulation and clarify the condition for the formulation that 
its solution satisfies the virial relation Madm = Mk. In Sec. IV, we derive a relation for the difference, SMadm — ^SJ, 
and clarify the conditions on the formulation for which a sequence of solutions satisfies the first law. In Sec. V, we 
propose formulations whose solutions and sequences of solutions approximately satisfy the virial relation and the first 
law. Section VI is devoted to a summary. 

Throughout this paper, we use geometrical units with G = 1 = c. Spacetime indices are Greek, spatial indices 
Latin, and the metric signature is — h H — h. Readers familiar with abstract indices can regard indices early in the 
alphabet as abstract, while i, j, k, - ■ ■ are concrete, associated with a chart {x*}. If 5* is a 2-surface in a 3-space T, and 
Cabc if' the volume form on S associated with a 3-metric 'yab, we write dSa = CabcdS^'^; for S a surface of constant r, 

dSa = VaTs/lCpX. 



II. FORMULATION 
A. 3-|-l formalism 

Let St be a family of spacelike hypersurfaces, labeled by a time function t. Let be a vector field transverse to 

for which f"Vai = 1, and denote by a and /?" a nonvanishing lapse and a shift vector, respectively, with 

t« = an" + f3", f3"na = 0. (2.1) 

Then, in a chart {t, x^}, we have = dt, and the metric g^fj = — n-an^ has the form 

ds^ = -a^dt^ + -iij{dx' + I3'dt){dx^ + P^dt). (2.2) 
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With a spatial covariant derivative Da compatible with the spatial metric 705, the extrinsic curvature of is given 

by 

Kab = -^£n7a6 = ^(-^«7a6 + D^Pb + DbPa), (2.3) 

where 'fab and dt^ab are the puUbacks to St of 7^/3 and £t7a/3- 

In the canonical formulation of general relativity [19], {706, tt"*", a and /3°} are regarded as independent gravitational 
field variables, where tt"'' is defined by 

tt"^ := -{K'''' - -i''^K)y/^. (2.4) 

A perfect fluid is described by a stress-energy tensor 

T"^ = (e + p)u"u^ + pg''^, (2.5) 

where p. and e arc the fluid four- velocity, pressure, and energy density, respectively. The pressure and the energy 
density are assumed to satisfy an equation of state of the form 

p = p{p,s), e = e(p, s), (2.6) 

where p is the baryon mass density and s the entropy per unit baryon mass. 

In calculating the variation of the Lagrangian following Routh procedure, a perfect-fluid spacetime is specified by 
the canonical variables, the lapse and the shift, that together describe the metric, and by Lagrangian variables for 
the fiuid, (5(A) := [7ab(A), 7r''''(A), a(A), /3"(A), ^"(A), p(A), s(A)]. The difference between two nearby solutions can be 
treated in either of two ways. Changes in the metric variables will be written as Eulerian variations, denoted by S; the 
Eulerian change is the difference between corresponding quantities in the two solutions at a fixed point in spacetime. 
Changes in fiuid variables will be written as Lagrangian variations. Introducing a Lagrangian displacement vector 
field one defines the Lagrangian change in any fluid variable as the change with respect to a frame dragged by 
Formally The Lagrangian change AQ in a quantity Q is then related to the Eulerian change 6Q by 

AQ = 5Q + £^Q. (2.7) 

The description of fluid perturbation in terms of a Lagrangian displacement ^" has a gauge freedom associated with 
a class of trivial displacements that yield no Eulerian change in the fluid variables. We use this freedom to choose a 
gauge in which ^* := ^"Vat = 0, following [20,21]. 
The Einstein-Hilbert action 



-I- 



Cd^x, (2.8) 

with the Lagrangian density 



takes, in terms of Hamiltonian metric variables, the form 

IGttC = 7r"*'5t7„6 - aHc - PaC^ + Da{-2D''a^ - 2(5\''b + /^^tt) - 9*77 - IStt ey^, (2.10) 
where is the Ricci scalar, 

WG:=-2G«Vn/3x/7, := -2G« V"/? %/7 , (2-11) 

and 0°^^ is the Einstein tensor. 

The variation in the Lagrangian density is given by 
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-SaH - 6P^Ca + dn''' \ dtjab - D^Pb - D^Pa - ^ (^TTab - llabTT ] \ - 5^ab{G'''' " 87r5"'')av^ 



-CaV0T«^^/^+£'„e«V7- Y^5t(^7r«Sa6) +a^(j„^^/7) , (2-12) 
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where T is the temperature, h is the enthalpy defined by h := (e +p)/p, and 

G"^= G"'^7„"7/ and S"^ ■.= T''^^^^^ . 



With definitions 



we set 



pjf:=T"Vn0 and f := -T'^'^ia^np 



7 2 J 



n := -2{G"I^ - 87rT«^)no,n^ ^ = Hq + l&ir Ph = 
C« := -2(G«'5 - 87rT«^)7a°n^ \/7 = - IGtt j«V7 = -2(£'67r«^ + 87rj«V7), 
where R is the Ricci scalar with respect to ^ab- The density is the surface term of the Lagrangian density, 



(2.13) 

(2.14) 

(2.15) 
(2.16) 



6" 
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- 25(£)«a V7) + {0''lhc5iT^' + 7r5/3" - 27r%(5/3'')} 



+ _ 7''''7-'')(ai?6,57cd - D^a <57cd) 



(2.17) 



where q"^ (g"'^ + u°'uP)-^^^i^ . 

Independently varying the metric variables, {5a,5(3°',5^ab,5i^°'^}, gives the field equations. 



W = 0, Ca = 0, and G"" - SttS"* = 



and the relation, 



dtlab - DaPb - Df,/3a - 



2a ( 



6 - -^abT^ 



0. 



(2.18) 



(2.19) 



Equation (2.19) is consistent with the definition of tt"'' [cf. Eq. (2.4)]. 

When the field equations are satisfied, the Bianchi identity implies W pT"-^ = 0. The variation of the action with 
respect to the (spatial) Lagrangian displacement vector is the spatial projection of this relation, the relativistic Euler 
equation. 



l^y^T^^ = 0. 

For an isentropic fluid, conservation of baryon mass and entropy are given by 

£u{p\/—g) = and £„s = 0. 



(2.20) 



(2.21) 



Equations (2.21) and (2.20) together imply V/^T"'^ = 0. 

It is often convenient to rewrite the above set of basic equations in terms of the conformally related spatial metric 
7o5 and the tracefree part of the extrinsic curvature Aab, defined by 



7ab ■= V' lab, 



Aab —Ip * (^Kab ~ ^labK^ , 



(2.22) 

(2.23) 



V, 



where is a conformal factor and K := K^b^y ■ Here, we may impose the condition, 7 := dct(7ab) = det(?7ah 
where r]ab is a flat 3-metric. In the following, indices of variables with a tilde, such as Aab, A°-^, (3a, and /?"(= /?°), 
are raised and lowered by 706 and 7"'', respectively. 

The Hamiltonian constraint W = and the momentum constraint Co = are written in terms of Kab'- 



R-KabK'^^ + K^ = 167tph, 
DbK\ - DaK = 8nja. 



(2.24) 
(2.25) 



With the conformal transformation (2.22) and (2.23), Eqs. (2.24) and (2.25) are rewritten in the form 
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(2.26) 
(2.27) 



Here, 7?, Da, and A arc the Ricci scalar, the covariant derivative, and the Laplacian with respect to 7abi respectively. 
The evolution equations for ^ab and Kab are 



dtlab= -2aKab + DaPb + DbPa, 

dtKab= aRab - DaDba + a{KKab - 2KacKi,' ) + {Db^Kca + {Daf3')K,b + {D^KabW 

Sab + -^lab [pH — •S'c^^ 



87ra 



(2.28) 
(2.29) 



where Rab is the Ricci tensor with respect to ^ab- 

Contracting 7'^'' with Eqs. (2.28) and (2.29) and using Eq. (2.24), one obtains 

dtK= aKabK"^ -Aa + 4wa{pH + SJ") + p'^daK, 



(2.30) 
(2.31) 



where A = DaD'^. To write the evolution equation of K in the form of Eq. (2.31), we use the Hamiltonian constraint 

equation (2.24). 

In the following, we choose the maximal time slicing condition K = = dtK. With this condition, Eq. (2.31) 
reduces to an elliptic equation for a, 



(2.32) 



where we keep spatial indices abstract until fixing the spatial gauge condition. Using Eq. (2.26), this equation is 
rewritten 



A(aV) = 2nai;^ipH + 25„") + ^a^^i^fei"*" + ^R . 



(2.33) 



Using Eqs. (2.28), (2.29), (2.30), and (2.31), the evolution equations for 7^6 and Aab are 



dtlab = -2aAab + DaPb + DbPa - -^labDcP" 



dfAab = ip~ 



a 



yRab - -Y^abR j - [DaDba - —7, 



abDcD'^a) 



-2aAacA^^ + bafi^A^b + b^P^A^ 
-87ra(^-45„6-^7a65/)- 

Now, Rab is split, 

Rab = Rab + Rtb' 

where Rab is the Ricci tensor with respect to jab and 



^b^fi^Aab+fi^b^Aab 



^ab 



2~ 6~~ 2 - ~ 

-DaDbiJ - -%b^^ + -n^Dat/jDbip - —.^bDc'^I^D'^^j:. 



1p' 



Rab is then written in the form 

Rab = l 



(0) (0) 



(0) 



(0) 



-T DcDd lab- Db {lacF^)- Da {ibcF") 

(0) (0) . (0) (0) . . ■ 

-{Dclbd) DaT''-{Dc lad) DbT'' + 2F'^C,,ab ' ^C^Cl^ 



' (0) - (0) (0) 

- A hab- Db ilacF^)- Da (ibcF^) 



, pNL 
^ab ' 



(2.34) 



(2.35) 
(2.36) 
(2.37) 



(2.38) 
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(0) (0) (0) (0) (0) 

where F" :=Db 7"^, -Da is the covariant derivative associated with rjab, and A= rj""' DcDd- Rat collection of 

the nonlinear terms in hab and defined by 
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■ . , (0) (0) . (0) (0) ~ , (0) . (0) . ■ 

DcDd hab + (Dc hbd) Da + {Dc had) Db 



+ F'^C'cab — C'^b^ad- 



Here, hab and /"^ are introduced, respectively, by 

Kb~7ab-Vab and /«^= ^^'^ - r,«^ 

C^(, and Cc,ab are defined by 

^-cd ^(0) _ (0) _ (0) . 



<^a6 ^bd+ Db had- Dd hab J and Cd,ab ■■= IcdCab, 



(2.39) 



(2.40) 



(2.41) 



(0) 



and C^^ =Dd Wj/v)/V^/'^ = 0' '^^'^en 7 = ?y. 



B. Basic equations for quasiequilibria 

Compact binary systems in quasicquilibrium circular orbits evolve toward merger due to gravitational radiation 
reaction. Since the emission time scale of gravitational waves is longer than the orbital period even just before the 
merger, we may expect that the fluid and field variables near the support of fiuid source (or inside the light cylinder) 
in the frame rotating with the same angular velocity as the orbital motion are approximately unchanged along a 
direction of a helical vector 

A:"=t" + r20", (2.42) 

where is a spatial vector field which generates a family of closed circular curves on St, and O denotes the orbital 
angular velocity. 

First, wc derive hydrodynamic equations to describe binary neutron stars in quasicquilibrium circular orbits. The 
baryon mass conservation law Eq. (2.21) and the Euler equation Eq. (2.20) are written 

,t 



£k+v{pu'y/^) = 0, (2.43) 

■ya"U+v{hUa) + Da (^^^ = 0. (2.44) 

Then, we impose the conditions that the Lie derivatives along fc" vanish 

£k{pu'V^) = 0, (2.45) 

Ja"£k{hUa)=0. (2.46) 

Here, a spatial velocity vector v°' is introduced by 

u" = u*(fc" +?;"). (2.47) 

From conditions (2.45) and (2.46), the relation £fc(.7aV^) ~ also follows. In the above, we assumed isentropic fiow, 
which leads to the local first law of the thermodynamics, 

-VaP = Vah. (2.48) 
P 

Equation (2.48) also implies that a one-parameter equation of state may be chosen. We thus have four independent 
variables for the fluid, three for the fluid velocity and one thermodynamic variable, governed by four equations (2.43) 
and (2.44). 

For corotational flow = u*k" (that is v" = 0) or irrotational flow hua = Va^, where $ is a velocity potential, 

one can obtain a first integral of the Euler equation that is useful for computing quasicquilibrium configurations. For 
corotational flow, the velocity field becomes trivial and the first integral is the statement that the injection energy is 
constant in the fluid: 
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h 

— = constant. 

For irrotational flow, one thermodynamic variable and a velocity potential are governed by two equations [22] , 



D. 



h ^ 
h_ 
7^ 



constant, 



0, 



(2.49) 

(2.50) 
(2.51) 



derived from Eqs. (2.43) and (2.44), respectively. Here, the Lie derivative with respect to the spatial vector v"', £„, 
is defined on St with a relation 



V" = Ta\ ^ -1^ 



(2.52) 



where is a rotational shift vector defined by = /3" + 00" . 

Note that the symmetry of the spacetime with respect to the helical vector fc" has not been imposed yet. Thus, the 
Lie derivatives of the fluid quantities along fc" may not vanish, e.g., £kP 7^ 0, £fe?ia 7^ 0, £kpH 7^ 0, and £kSab 7^ 0. 
The values of these quantities depend on the formulation for gravitational fields that we choose below; their magnitude 
measures the deviation from the helical symmetry. 

We turn to the formulation for the gravitational fields of binary systems in quasicircular orbits. Here, we do not 
assume a global helical symmetry for the whole spacetime. First, we define u"'' :— dtj"''' and regard it as an input 
quantity: It is determined when we impose a certain condition between 7"^ on two spatial hypcrsurfaccs of infinitesimal 
time difference, following the concept of a thin sandwich formalism proposed by York [23]. In this section, we continue 
calculation without fixing the condition for except for a requirement 

u"^ = dtT^ = 0{r-'^), (2.53) 

in a far zone (r 27rr2~^). This condition guarantees the asymptotic flatness of the system on a slice St, but breaks 
the helical symmetry for 7"^ in the far zone. The lower component of the time derivative is defined by 

,,cd 



Uab : = 



-laclbdU 



We also define 



Vab ■■= dtAab, 



(2.54) 



(2.55) 



where Aab ■= ip^Aab- Later, we will impose a condition on Vab- 

Bonazzola et al. [15] propose a "gravito-inelastic approximation" in which they set u"'^ = 0, while dtA°'^ is free. On 
the other hand, Schafer and Gopakumar [16] propose a minimal no-radiation approximation, in which the transverse- 
tracefree part of dtir"'^ is restricted. We determine the conditions in the more rigid way: Our conditions for u"'^ and 
Vab are determined from the requirement that the quasiequilibrium solution and its sequence satisfy the virial relation 
and first law at least approximately. This subject will be discussed in Sees. HI and IV. 

Equation (2.34) is regarded as the equation for determining Aab, namely, 



2aAab = D^Pb + DbPa - T;%bDcf3' 



Uab- 



(2.56) 



Then, substituting Eq. (2.56) into Eq. (2.27), we obtain 



A/3„ + ^DaDbP" + RabP' + In ( ^ 

a 



DbPa + Da(3b - -labDcP' 



^b,(a-^ilj^f'Uab) = 167rai„. (2.57) 



This elliptic equation determines Pa- 

Regarding Vab as an input quantity, the evolution equation for Aab may be rewritten as an elliptic equation for hab, 



(0) 

A hab ■ 



Rtb - \ Da {F^%b) - \ Db {F'%a) + R^^ - ^labR ' \ (DaDba - ^7a6Aa) 



Ai^'^AacAC + — ( Dal3'A,b + Dbf3'Aac 



a 



^Del3'Aab 



- P'DcAal 



(2.58) 
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Together with the above equations for hat and 0a, the elhptic equations (2.26) and (2.33) are solved for ip and % := aip, 
respectively. 

To summarize, the Einstein equations are rewritten as four elliptic equations (2.26), (2.33), (2.57), and (2.58) for 

ip, X: and hab, with m"^ = 9(7°'' and Vat = dtAat which are regarded as input quantities that satisfy a certain 
ansatz. The asymptotic behavior of u"'^ and Vab is chosen to preclude standing waves in the far zone, r ^ 2Tr^l~^. 
The simultaneous equations for the metric and the fluid on a slice Sj is then similar to the equations for the initial 
value problem since the equations for the metric variables are elliptic. 

To solve the equations for hab, we need to fix the spatial gauge. The simplest choice is a transverse gauge [4] (or a 
generalized Dirac gauge in the terminology of [15]), satisfying 

(0) 

F" =DbT^O- (2.59) 

With this choice, the equations for hij and R are significantly simplified, resulting in that the operator for the linear 
terms of hij in its elliptic equation becomes the flat Laplacian. Furthermore, the behavior of the source terms of the 
elliptic equations for ip and x for r » 27rr2~^ have suitable asymptotic behavior, because = 0{r~^) and, hence, 
R = 0{r^*). Thus, the source terms of elliptic equations for %b and \ arc 0(r~^) in the present gauge. This implies 
that we can numerically solve these equations without serious difficulties. 

As a result of the present choice of gauge conditions, the asymptotic behavior as r — > oo of the geometric variables 
is [24] 

^|J = l + ^^+0{r-% (2.60) 
2r 

a = l-M^+0{r-^), (2.61) 



r 

= -^(^eZkif' + SZijf'f^r'' - Zur'^ + SZ^i'r'^ + 0{r-^), (2.62) 

hij = 0{r-^), dkhij = 0{r-^), dkdihij = 0{r-^), (2.63) 
1 r 

6Zij - 26ijZkk - QZiir'-f^ - QZjif^f' + ASijZuf^f'' 



A --^ 



+ [bZkifH' - Zu){5ij - 3f'F) - UiZf^^f'P + Zf^rH') + U,, 



+ 0(r-4), (2.64) 



where Madm and Mk denote the ADM mass and Komar mass (see Sec. Ill), = /r, and Z^j and Z^^ are 
time-dependent symmetric and antisymmetric moments, respectively. Equation (2.64) implies that the total linear 
momentum of the system is implicitly assumed to be zero, as 

/ AldSi= [ i)^A/dSi=0. (2.65) 

J oo J oo 

Here 

/ := lim / , 

Joo JSr 

with Sr a sphere of constant r. The antisymmetric moment can be related to the angular momentum of the system 
as 

= -Jj<^jki, (2.66) 
where ejki is the completely antisymmetric symbol. In the Newtonian limit, 

Zki = I piv'x' + v'x'')Sx = (2.67) 

Z^^ = J p{v''x^ - v^x'')d^x, (2.68) 

where lij is the quadrupole moment. Uij in Eq. (2.64) is a symmetric tracefree moment determined by the asymptotic 
behavior of u,- 
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A solution to the simultaneous equations derived here satisfies the constraint equations of the Einstein equation. 
In this sense, it can be referred to as a fully general relativistic solution and can be used as an initial condition of the 
(3+1) numerical simulation. However, since we do not assume the helical symmetric relation £kjab 7^ 0, the solution 
does not satisfy the equation for quasiequilibrium exactly. As a result, £kij, £-kce, and £/c/3° would be slightly different 
from zero in general. Deviation of these quantities from zero can measure the violation of the helical symmetry. The 
magnitude of the deviation depends on our choice of u"^ and Vab- 

III. RELATIONS BETWEEN Madm AND Mk 

In this section, we derive the conditions needed for equality of the ADM mass Madm and the Komar mass Mk- 
This equality is closely related to the virial relations as discussed in Appendix A. 

1. Sufficient conditions for equality of Mk and Madm 

The Komar mass [25] is constructed from a vector ^" that approaches a timelike Killing vector of a flat asymptotic 
metric at spatial infinity. As presented in [26], = — Q!^V"t, and 

= ^ y {v^c - v^c^wdSff. (3.1) 

With the metric written in the form (2.2), one uses the relations = aW at and n^W pn" = a~^^^^a(^-, to obtain 

V/3(V''C" - V«C'')n„ = aV/3{7^[V«(an^) - VT(an")]Vat} = 2£)„D«a. (3.2) 
Using Eq. (3.2), we have 

i- j 2DaD''adV =^jD''a dSa, (3.3) 
and thus, we reach the familiar form [27] 

Mk = ^ f D'^a dSa. (3.4) 

For a stationary spacctime, with C" the asymptotically timelike Killing vector, Beig [28] and Ashtekar and Ashtekar 
[29] prove the equality Mk — Madm- Wc obtain here more general asymptotic conditions sufficient for equality in 
the following way (patterned in part on Beig's work). Suppose that the metric has the form (2.2), with^ 

lij =Vij (3.5) 

hij = 0{r-^), Dkhij = o(r-i), DkDihij = o(r-3/2), (3.6) 
and suppose that the lapse, shift, and extrinsic curvature satisfy 

a = l-- + o(r-i), dra = ^ + o{r-^), d^a = -2^ + o{r-^y, (3.7) 
dQKrr = o{r-^), Kij=o{r-'"% dkKij = o{r-^+'), (3.9) 



L 



where / dCl denotes a surface integral. Then Mk = Madm- 



^The definition hij here is slightly different from hij in Eq. (2.40) in Sec. II B. 
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To prove this claim, it is useful to introduce S^Gij, the part of ^Gij linear in hij. The idea is to show that if S^Gij 
is the asymptotically dominant part of ^G^ , then 

/ ^Gix^ dSi = -SttMadm. (3.10) 

J oo 

One then uses the field equation for ^Gij (the dynamical equation for Kij) to show that this integral can be written 
in the form (3.4) of the Komar mass, when Kij, £fjKij, and dtK^r fall off rapidly enough at spatial infinity. 
Formally, 

, 1/(0) (0) , (0) (0) (0) (0) (0) (0) (0) (0) N 

5 Gij = -(^ DiDk h]+ DjDk hf- A hij- DiDj hl+rjij A K - rjij DkDi h^^j, (3.11) 

where the index of hij is raised by the flat metric t?'-' . The asymptotic behavior of ^Gij is given by 

^Gij - S^Gij = o{r^). (3.12) 

This is because each term in ^Gij involves either DkDihij or DkhijDihmnl then terms quadratic or higher order in 
hij fall off as rapidly as either hmnDkDihij or Di~hijDihmn- 

(0) 

From the linearized Bianchi identity, D-'o°Gij = 0, we have 

j ^G'jX^dSi = J D\5^GijX^)^(fx = J 6^Gijrj'^,/^(fx, (3.13) 



where 



Then 



1 (0) (0) 1 (0) , 1 (0) (0) , (0) 

rj^^S'Gij = -- DiDj h'' + -Ahl = - Di {D'h^ - Dj h^'). (3.14) 
/ ^Gja;^' dSi = -\ I {Dj h'^- D'h'l)dSi =: -SttMadm- (3.15) 



Next, from Eqs. (2.29) and (2.24), Gij has, outside the matter, the form 
Gij = Rij — —JijR 

= -DiDja - l:J^J-Aa + -(OtK^j - \jijj''^dtKki) + 2KikK^ - KijK 

-7y (KuK'^' - i^^^ + ^ZnKij - ^ji^^'^'E^Kki. (3.16) 
Since the terms involving Kij are o(r~^), we have 

/ ^G'jX^ dSi = j DiDja f'f^r^dfl. (3.17) 

J oo J oo 

Now dkjij = 0{r~^) and dia = 0{r~^) imply 

Dj^a = dj^a + Oir'^). (3.18) 
Finally, from Eq. (3.4) and the asymptotic form (3.7) of a, we have 

Mk = — [ hidh, (3.19) 

and 

/ ^Gja;-' dSi = ( D^^a dCl = -SttMk, (3.20) 

J oo J CO 

whence 

Madm - Mk. (3.21) 

The asymptotic behavior of Eqs. (2.60)-(2.64) does not always satisfy the conditions (3.9), because dtKij is 0{r~^), 
not o(r~^). 
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2. Relation between Mk and Madm in the waveless approximation 



For quasiequilibrium binary solutions that satisfy Eqs. (2.54)-(2.58), with K = 0, a, slightly different relation holds 
between Mk and Madm- The asymptotic behavior given in Eqs. (2.63) and (2.60) implies 

Dihij = 0(r-2), DiDjhij = 0(r-3). (3.22) 

For the asymptotic behavior of Eqs. (2.60)-(2.64), together with Eq. (3.22), Eq. (3.15) still holds and 

/ ^G^x^dSi = -SttMk + / dtKrydh = -SttMk + j dt{l2Zrr - Urr)dh, (3.23) 



or 



dt\2 

Using 



Madm = Mk - ^ {\Zrr - \Urr ] ■ (3.24) 



/ 



f'f^dn=^5ij, (3.25) 



Equation (3.23) may be written in asymptotically Cartesian coordinates as 

/ ^G]x'dSi = -SttMk + / ^tXy fV^r^dJl = -SttMk + 47r^^, (3.26) 

whence 

Madm = Mk-^^. (3.27) 

We note that a similar expression has been derived for the maximal slicing condition in the first post Newtonian 
approximation [30] . The expression here is the fully general relativistic generalization. 
In the presence of a timelike Killing vector, we may set 

dt%b = 0, (3.28) 

at(^V'') = 0, (3.29) 

j^Zkk = 0, (3.30) 

and, hence, the virial relation Madm = Mk holds. This condition should be satisfied not only for axisymmetric 
equilibria but also for nonaxisymmetric ones such as general relativistic Dedekind solutions. 
The integral of Kij in the above expression may be computed further as 

I ^^^"^ =ljtL ^'^'^'''^ -Ijtl 

= -Jt I '^^^'""^WA I -"-^SU..^^., (3.31) 

which yields 

Madm = + j„ dV - ^ j {dtn"'' Dc Jab + tt"^ Dc daab)d''x 

" + 1 / " 1^ / [^*(^'''"') ^-^^ + daab]d^x. (3.32) 
In a gauge with K = 0, Eq. (3.32) can be written as 
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Madm = Mk + j^J x'^jadV - [dtAabx'^Dc r'' + Aabx''Dc dtT'']d^x 

■ j X'^ja dV - [Vab X" Do T'' + ^ab x" Do u'''>]d^X, (3.33) 



' dt 

(0) , (0) 

where a relation tt"'' x'^ Dc lab = -^ab x'^ Dc 1°'^ is used. Since we always impose a condition £k{ja\/l) = 0) the second 
integral can be discarded because of the relation, 



x'^jadV = / dt{x''jaVl)d''x = / £k{x''jaV^)d''x = 0. (3.34) 



d_ 
di 

Thus, Eq. (3.33) is written as 

Madm = Mk - / {vab x' Dc T'' + ^6 x' Do u'''>)d^x. (3.35) 

If we assume that the spacetime is everywhere helically symmetric, £kgaf3 = 0, then it is not asymptotically flat 
and, hence, the integrals in Eqs. (3.32) and (3.33), as well as Madm, are not defined. Instead, we require helical 
symmetry only in the near zone and the local distant zone. In this case, in the local zone, we should set £fe7"^ = 
and £kAab = which can be written as 

= atr'' = -£o0r', (3.36) 

Vab = dtAab = -£n4,Aab, (3.37) 

while in the distant zone, w"^ and Vab ^ for r ^ oo. With these choices, dZkk/dt vanishes, and hence, the 
virial relation Madm = Mk is satisfied. The virial relation is satisfied for the simple case u"'' = Vab = 0; Equation 
(3.35) reminds us that it is satisfied in the conformal flatness approximation 

^ab^^ab^ i.e., U«^ = [21]. 

IV. RELATION FOR SMadm AND dJ 

A first law for binary systems with a helical Killing field k" has been formulated as relating a change in a conserved 
charge Q, associated with a family of helically symmetric spacetimes [21], to the changes in the vorticity, baryon mass, 
and entropy of the fluid as well as in the area of black holes. In [21], we prove that a relation 6Q = SMadm — ^SJ 
holds for asymptotically flat systems, where is the orbital angular velocity of a binary system in circular orbits. 

In inspiraling binary neutron stars, entropy, baryon mass, and vorticity are almost constant and, hence, energy 
and angular momentum are dissipated only by gravitational radiation. Thus, the relation dMABu/dt = fldJ/dt is 
satisfied. Then, we should require that, in a formalism for computing asymptotically flat binary equilibria, a first law 
SMadm = ^SJ is satisfied. In this section, wc present a heuristic way to derive a relation held between the variations 
of the ADM mass and the angular momentum without assuming any symmetry for perfect-fluid spacetimes, but 
requiring the fleld equations to be satisfled. Then, we identify sources for the violation of the flrst law. In the 
following calculation, no gauge condition is specified; and surface terms associated with black hole horizons are not 
included. 

In contrast to our earlier paper [21], no helical symmetry is assumed. Instead, we assume only that the field 
equations derived in Sec. II are satisfied. The field equations are used to relate the Komar mass to the Lagrangian 
density and the other terms in the manner 

= ^ ^ D'^adSa = ^J D'^D.adV 

= 2 J i^pu'^huc^v" -C- A.^^^a^TT"'' - ^£)„(27r»V6 - 7r/3»)| d^'x, (4.1) 

where the definition of v°', slightly different from the previous section, is = u*{t°'+v°') with v"na = Q- To compute 
the last equality, we first use a part of the Einstein equation, 

(C?«^ - STrT"'^) (^7«^ + Ua^np + ^/^.n^) 
= 7a6St7r"^ + ID^DaOL^ + 167r {A^^ - \/^ - 167rpw*v/^/iw„t;« + Da{2TT"^(3b - tt/?") = 0. (4.2) 
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We then subtract the trace of the Einstein equation (G"'^ - 8iTT"^)ga0 = -R- SwT^" = to relate Mr to the 
Lagrangian density (2.9). 

The angular momentum is defined by 



where the momentum constraint C = is used in the last equality. 

The variations of Mk and J are computed following [21]. First, we take the variation of Eq. (4.1) for Mk, 

and substitute the following relation into the first term, 



(4.3) 



(4.4) 



(4.5) 



where we used the facts that a choice of ^* = implies Au" = —At" = £t^" and that £t^"Va^ = yields 

pu^ ^—ghua£,tSf = jay^dtS,"'. To the second term, the variation in the Lagrangian density Eq. (2.12) is applied. 

When all components of the Einstein equation and the Bianchi identity, Eqs. (2.18)-(2.21), are satisfied, 6Mk 
becomes 



(4.6) 



The contribution of the surface term of the above equation is given by 



-2 jo. 



e"^/7+Y^'^(27r"V6-7r/?») 



d^x = -2 / ^"dSa 

J oo 



-3- / D-6adSa I (7"'^7'" " 7"'7''')£'b<^7cdd5a = <5Mk - 2(5Madm. 

47r 87r 



(4.7) 



Combining Eqs. (4.6) and (4.7), the variation in the ADM mass Madm, instead of Komar mass Mk, is written 



(5Ma 



DM 



pTAsV^ +(— + hu^vA A{pu*y^) + v"A{hu„)pu'^ - ^dtUaVl) 



+ (^7r«''at7a6 - <57a 

167r 



yabdtn"'')^ d' 

The variation in the angular momentum 6 J, computed in a similar way from Eq. (4.3), is 
6J = j |A(i„<^«V7) - Y^^(7r"''£07a6)} d'^x 

= y I hua^rnpu'v^) + r^ihu„)pu'v^ + r£0(jaV7) 



(4.8) 



-^(<57r«^£^7a6 - SjabS^^T^'") ' U ( jaC + T^^^'^^ab 1 \ d'x, 



-ah 



167r 



167r 



where 



A(jar V7) = A(p«*\/^)/iu„0" + pu'^rA{huo) - iaV7£0r 

= A{pu'^)hUc.cj)'' + pu'^^'^Aihu^) + r£0(.?aV7) - £0(.?arV7), 



(4.9) 



(4.10) 



which results from relations ja</'"V7 = pu* \/—9hua4'°' , 6(j)°' = 0, and Ac/)" = —£0^". The last term in the integral of 
Eq. (4.9) vanishes since it becomes a surface integral with a combination (j)°'dSa = 0. 



13 



Finally, Eqs. (4.8) and (4.9) are combined to derive a relation 

h 



SMa_ 



DM 



I 



pTAsV^ - 



+ hua{v" - 00") \ A(pu*v^) + {v°' - f20«)A(/iWc«)/9w*y^ 



(4.11) 



where £t+Q4> = dt + £r20, which operates on the spatial quantities, is understood as a puUback of the Lie derivative 
along the helical vector field onto St. The first three terms of Eq. (4.11) are the same as those derived in our previous 
paper [21] except for the definition of the spatial velocity v" = u"/u* (in [21], we used the definition i;" = u°'/u* — Q(j)°'). 
The difference also changes the definition of the shift. 

As discussed in [21], for an isentropic fluid, conservation of baryon mass, entropy, and vorticity become 

£u(p\/^) = 0, £„s = 0, and £„a;c«/3 = 0. (4.12) 
These imply perturbed conservation laws, 

A{pu*y^) = 0, As = 0, and Aw^/s = 0, (4.13) 

that will be almost satisfied during binary inspiral before the merger. Here, the relativistic vorticity is given by 

= qa'^q/}^ [V^ {hus) - Vs {hu^)] = Va {hU/}) - V/3 {hUa) . (4.14) 

The third term in Eq. (4.11) vanishes for (i) corotating binaries, flows with v" = f2<^", and (ii) irrotational binaries, 
potential flows with hUa = Va$. For the latter case, the third term in Eq. (4.11) with A{hUa) = AVa^ = VaA$ 
becomes 

j {v°' - ^(?i")A(/lu„)pu*^/^rf^a; = j (v" - n^")VaA^pu*s/^ d^x 
= j [£fc(py^) - £„(p^)]A$d3a;, 



(4.15) 



where we assume p = for the distant zone to derive the last line. Then, together with Eqs. (4.12) and (4.13), Eq. 
(4.11) is rewritten as 



5M, 



ADM 



(4.16) 



The form of the "first law" described in Eq. (4.11) or Eq. (4.16) is derived without relying on the helical symmetry 
of the spacctimc and the fluid. Choosing the maximal slicing condition if = = tt, we may rewrite Eq. (4.16) using 



tracefree part of the conjugate momentum t:""^ = n"'^ — -^""^tt and the conformal metric ^ab ~ ^7a6 as 



i5Madm - = 



£fe(pV^)A$ - ezkijaVi) + -7r{5{^ ^ )^kiab - 5^ab£k{r*n} 



d^'x 



-j 



£fc(px/^)A$ - ^£fe(ia^/7) + j^i^Aab^kT" - ST'^kAab) 



d-'x. 



(4.17) 



Here, the gauge choice tt — Q implies 5-k = Q and £fc7r = 0. 

If one requires £fc7a6 = and £fc7r"'' (or £^7"*' = and £fe^ab = 0) in a gauge n = = K, together 
with conditions for the fluid variables £k{p\/—g) = and £fe(jaV7) — 0' Eq. (4.16) leads to the first law relation 
(5AfADM = i^SJ. However, these assumptions are equivalent to imposing helical symmetry on the whole spacetime 
and, hence, preclude asymptotic flatness; in other words, Madm and J arc ill defined. 

In a realistic system, the radiation reaction violates the Killing symmetry. In its presence, Sjab and Jtt"'' are 
determined by the radiation reaction, and these terms may be proportional to the violation of the helical symmetry 
near the source. Namely, we expect that the following relations hold: 
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{£k7ab)St = (4.18) 
(£fc7r"'')(5t = Sir"'', (4.19) 



or in the gauge K = 0, 

{£kr'')St = Sr\ (4.20) 
{£kAab)St = SAab, (4.21) 

where 6t is a radiation reaction time scale. In this case, the right-hand side of Eq. (4.16) vanishes with the symmetry 
for the fluid variables. This indicates that even with the slight violation of the helical symmetry due to the radiation 
reaction, a relation SMadm = ^SJ may be well satisfied. 

Finally, as shown in [21], SMadm = ^SJ is exact in the conformal flatness approximation (IWM formalism). In 
this case, one needs to replace the Lagrangian density (2.9) by one that reproduces the field equations of the IWM 
formalism. One can derive such a Lagrangian density by substituting tt = and 7ah — rjab into Eq. (2.9). Then, 
assuming helical symmetry for the fluid and from the fact Sjab = 0, the first law is shown to be satisfied. (See [21] 
for a description of the artificiality of this choice in a helically symmetric IWM framework.) 



V. CANDIDATE FORMULATIONS FOR QUASIEQUILIBRIA 

The condition u"^ = 0{r~^) is not compatible with helical symmetry in the whole spacetime. Thus, we propose to 
impose 



-£0^7"" forr<ro, 
for r > ro, 



(5.1) 



where ro is an arbitrary radius. With this condition, the type of the field equation for hij changes from Helmholtz-type 
to elliptic for r ~ tq. To make the equation be almost elliptic for numerical computation, it may be desirable to take 
ro within the light cylinder radius as ro ^ 2tt/CI. On the other hand, we can impose helical symmetry on Aab without 
serious difficulty. In this case, helical symmetry is exact in the near zone and, as a result, the violation of the first 
law is given by 



^Madm -n5J = -^ f {6Aab)£n<i>r''d''x. (5.2) 

iDTT ./r>ro 



Since (<5Aaf,)£o07°'' falls off as 0(r^^) and the integral is done only in the distant zone, the magnitude of the integral 
would be very small. Thus, even with the modified formulation, the first law would be satisfied approximately. 
Furthermore, the virial relation is satisfied in this formulation. 
The condition for dtAab may be changed to 

va, = i -^"^^'^^ ' i (5.3) 
10 for r > ro. 



Then, 



(5Madm - ^SJ 

I I-ITT f 

r>ro 



— / 



(<5Aab)£n^7'" - (£n0^a6)'57"'' 



d^x. (5.4) 



Even in this case, the magnitude of the violation of the first law would be small, and the virial relation holds. The 
merit in this approach is that the right-hand side of the elliptic equation for hab falls off as 0{r~^). As a result, it is 
numerically easier to integrate the equation. 

We also note that instead of using the step function, we may write 

u<''> = -£j^r)mr\ (5.5) 

Vab = -£j{r)n0Aab, (5.6) 

where /(r) is a smooth function that satisfies the condition 
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1 for r <C ro, 
for r » ro- 



This choice is equivalent to taking a Kihing vector of the form 



dt 



+ /(r)0 



(5.7) 



(5.8) 



This Kilhng vector is hehcal in the near zone and purely timelike for r oo. 

Finally wc comment on other possible formulations. In the formulation with = = Vab, the virial relation 
is satisfied for a quasiequilibrium binary. However, the first law along quasiequilibrium sequences is not satisfied in 
general. The violation of the first law is written as 



smaum -nsj = — 

Ion 



167r 



(5.9) 



VI. SUMMARY 



Two relations, the virial relation Madm = Mk and the first law 6Mabm = ^SJ, are regarded as guiding principles 
to develop a formalism for computing binary compact objects in quasiequilibrium circular orbits in general relativ- 
ity. Deriving the explicit equations for Madm — Mk and SMabm — ^SJ on the assumption that the spacetime is 
asymptotically flat, it is shown that a solution and a sequence of the solutions computed in some formulations satisfy 
these two conditions at least approximately. We propose a formulation in which the full Einstein equation is solved 
with the maximal slicing and in a transverse gauge for the conformal three-metric. In the proposed formulation, the 
solution in the near zone is helically symmetric, but in the distant zone, it is asymptotically waveless. 

So far, quasiequilibria of binary neutron stars have been computed using the conformal flatness approximation for 
the three-metric [9,11]. In this formulation, only five components of the Einstein equation are satisfied, and thus, 
the obtained numerical solutions for quasiequilibria involve a systematic error. Specifically, in a real solution of the 
quasiequilibrium circular orbit, the conformal nonflat part of the three-metric will be of order (M/a)^, which can be 
~ 0.1 near the neutron stars for close circular orbits of a < lOM (e.g., [4]). This implies that to compute an accurate 
quasiequilibrium in circular orbits of error within, say, 1%, it will be necessary to take into account the conformal 
nonflat part of the three-metric. In the new formulations described here, such term is computed, and thus, more 
accurate solutions of quasiequilibria will be obtained. Currently, we are working in computation of binary neutron 
stars in quasiequilibrium circular orbits using these formulations. In a subsequent paper [31], we will present the 
numerical results. Such a numerical solution will be also used as an appropriate initial condition for simulations of 
binary neutron star mergers [5]. 

In this paper, we restrict our attention to the system in which no black hole exist. In the presence of black holes, 
we should carefully treat the surface terms at event horizons. The surface terms would modify the equations for the 
virial relation and first law [21,33]. The formulation for computation of quasiequilibrium black hole binaries are left 
for the future [32]. 
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APPENDIX A: EVOLUTION EQUATION FOR THE SCALAR MOMENT AND VIRIAL RELATION 

In this section, we derive the virial relation by direct integration of the Euler equation. Thus, the virial relation we 
consider here is associated with an evolution equation for the scalar moment as in the Newtonian case. In the end, 
we confirm that the virial relation derived is equivalent to Mk = Madm- 
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In the following, we often refer to -M^ as a "Komar like mass" that is defined by the asymptotic behavior of a 
function x '■= Q^V' at r ^ oo, 



X-l-^ + 0(r-^). 



(Al) 



For simplicity, in the following calculation, wo adopt a gauge in which if = 0, — 0, and j — r/, and we carry out 
the calculations in Cartesian coordinates. Wc often need to evaluate surface integrals at r ^ oo. In the evaluation, 
we assume Eqs. (2.60)-(2.64) as well as (Al) the asymptotic behaviors of geometric variables. As a consequence, all 
the volume integrals that appear below are well defined, and furthermore, the surface integrals derived during the 
calculation can be safely discarded. 

From the asymptotic behavior as r ^ oo, we can define and Madm using the surface integrals 



Mv = 



2tt 



Madm 



Using the Gauss's law, they can be rewritten in other forms 



2tt 



Maum 



1 

'2^ 



{xAij + f^diijdjx)d''x. 



The difference between Madm and is written in the form 
- Madm = 



d^x. 



Here, using Eq. (2.34), we can derive an identity, 

j ai,^A'jA\d'x = i y +7ifc7^"a„/3'= +y"/3'=afc7m - ¥''uik)d'x 

[ - {dii^'A'^) - ^'A'ki'ji}p' + l^'^Aiju 

[Snji^'^P' -^A^ju'^jd^x, 

where ^^j denotes the Christoffel symbol with respect to 7ij , and Aij = ijj^Aij . Thus, we obtain 
Mx - Madm = j 



I 



d^x 
d^x 



iDTT OTT TT 



d-^x, 



where we use the relation 



= 



kl 



+ 3aP = jk{v''+p'')+3aP. 



(A2) 
(A3) 

(A4) 
(A5) 

(A6) 



(A7) 



(A8) 



(A9) 



(Note that v'^ here is defined by v'^ = u''/u*.) In stationary spacetimes, the relation Madm = M^^ [28,29] and 
Uij = dt'jij = should hold. Thus, we get the virial relation as 



i;^ji{2v' - (3') + 6ai;^P + ^xi^R + -f■■'{^^^)^JX 

OTT TT 



d^x = 0. 



(AlO) 



In quasiequilibrium binaries, Uij 7^ in general. Thus, the virial relation, Madm = M^, is written as 
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IGtt 



1 x^^+-f^(aiV)ajx 



d^a; = 0. 



(All) 



As in the Newtonian case, we can derive the general relativistic virial relation from the evolution equation for the 
scalar moment. First, we write the general relativistic Euler equation ^"]^ ixT% = in the form 

StOfe^^) + ^ a^(av;6p) + p^^^d^x - [PH + 2Si')x4'^du^ - i^^'jAfi' + \xi>s,,d^f^ = o. (A12) 

Equation (A12) is a fully general relativistic expression, and no simplification is done. Taking an inner product with 
a;*^, we have 



In the following, we carry out the integral for each term separately. 
(1) First term: 



(fx = 0. 



dt 



x'^jki/j'^d^x. 



(A13) 



(A14) 



In the Newtonian limit, jfeV^ ~^ P^*^ = pdx'^/dt, and thus, this term leads to half of the second time-derivative of the 
scalar moment, i.e., /fefe/2 . 

(2) Second and third terms: By integration by parts, we immediately find 



h:^ J x''d,{jkv'ind:'x = - j Jkv^^j'd-'x, 

h:= j x''dk{ai)^P)d^x = j atjj^Pd^x. 

In the Newtonian limit, —I2 and —Is are the terms associated with kinetic energy and internal energy. 
(3) Fourth and fifth terms: Using Eqs. (2.26) and (2.33), we can rewrite the combination of them as 



R 



1 



PHrdkX - (PH + 2Si')xrdu4' = —dkixi^) - — 



167r 



27r 



(AV)afeX+(Ax)afeV 



167r 



Taking into account an identity, 

J [{x''dk^)Ax + [x^dkx)^^ 
where we discard the vanishing surface integral terms, we find 



/ 



l'^{dix)dji> + x\dix){dji>)dkl'^ d^x, 



x^'M'dkX - {PH + 25,')x^"afe^]d3 



1 

16tt 



(4) Sixth term: 



h:-- 



Rx'^dkix^) - ^[l'\diX)dji} + x\dkf^){dix)dji}] - i'^^A/A/x'^dk 
-- - J ilj''jix''dkP'd^x 



d'^x. 



-^^A^^x'^OkdiP' + diP') + -x\dkP')i>''Aijdif^ 



1 

"8n 



Stt/ 



{d,p')x''dk{i^''A\) + 2(5,/3V'^i' + -x''{dkf3')i^'A,,dif^ 
{dip')x''dk{i^'A\) - 2p'di{i^'Ai') + ^x''{dkP')^^Aijdif^ 



d'^x 
d^x. 



(A15) 
(A16) 



(A17) 



(A18) 



(A19) 



(A20) 
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Here, let us evaluate the first term. Using Eq. (2.34), 



{2aA\-j^'(3''dk%+7''uij)x''dk{ij''A\) + i^Jp'')i:''A\x"'^m{J''J^k) 



1 

"Stt 



(fx-lL 



(A21) 



Thus, 



p — _ 

167r 



167r/ 
1 

~ ~167r 



As a result, 
h-- 



{2aA\ - ^'=|3^^kl^J + l''u,,)dra{ll^''A\) + {^,|3''){A''^ml^k + AklOral^'W 



x'^d^x 



dm{i^'^A'iA\) + {%P''dk^^'+7''uij)dm{i^'A\) + {dj0''){A'^dmlik+Audml'')i^' 



x'^d^x 



d'^x. 



1 

167r 



j -^ai^'A\A\ - i^'^A\A\x"^d^ (^^) + ^\dki^'){dm{'^^ A^i) - i^'A\dmli,)x" 



+x\di,l3')i,''A,,dif^ - 2/3'(167rizV' - i'^'A^^dif^) 



+x\dki3')rA,,da'^ - P^Snji^p'^ - 2^p^A,,dif^) 



d'x. 



(5) Seventh term: Using Eq. (2.35), we rewrite it as 

I iOTT [V ■'a 

- 2aiifei/ + bi'fi^Akj + Dj~0^Au - ^Dk^'^Aij + (3^DkAij - dtA, 



where we use jijdkY'' = ln7 = 0. 

By straightforward calculations, we obtain 



Ir-.-- 



167r 



dkf'd^x = — I [- DiDjixi^) + m^)djx\ x^duf^d^x 



167r/ 



(A22) 



(A23) 



(A24) 



h~ j x^x\dkf^)Rijd''x = - / [di{xi^)x\dkf')T\^ + x''dk{xi')R + Rxi] dPx, (A25) 



(A26) 



+ m^)djx\x''dkf'd''x. 

Here, to derive the first equation, we use the spatial gauge condition F'^ = and relations in the present gauge as 

(A27) 
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The spatial gauge condition is also used in calculation for /y. 

To evaluate the remaining terms, we first rewrite the following equation using the definition of Aij as 

= p'^dkAij - A/^udkP' + AikdjP" - A'jP'^dkiu + A'jUik. 
Then, after a straightforward calculation, we get 



(A28) 



/8= J -2aAikA''j + bip'^Akj + DjP'^Aki - ^DkP'^Aij + p'^DkAtj - dtA^ 
"16^/ 



+ i,^x^A'^{dkl3^)di^ij + x^i}^Aiu{djl3^)difi - /3^x"V'i'i(9fc7i;)5n7'^' + x^^"" A'^^uM'- 



(A29) 



where we use an identity dtlmj)^ = Dk(3^ = ip ^dki'tp^P'^) that follows from the maximal slicing condition K = 0. 
Eventually, we find that I5 + h has the following simple form 



1 

167r 



Akiu""' +i;'^A/A,'x^d„ 



d-'x. 



(A30) 



By summation of /i ~ /g, we obtain the following simple relation: 



d'^x 



Madm - My d_ 
2 dt 



J x'^jkip^d^X-^ J {VijX^dnf' + AijX^dnU'^)d^X 



Here, since = dtj^^ and Vij = dfAij, the second integral term in the last line of Eq. (A31) is rewritten as 



1 d 

16^ di 



{^''A,,x''dr,f^)d^X. 



(A31) 



(A32) 



Using the momentum constraint, we further rewrite this term as 



"SiTdi 
d_ 
"dt 



d,{A\^^) + -^^Aijdkf' \d'x 



I 



dSiA\x''iP^ 



I 



x^juVd'^x - -Zkk 



Thus, a similar relation between Madm and M^ 



Madm = M^ 



dt 



(A33) 



(A34) 



is derived as is done for Madm and Mk in Sec. III. In this way, one can associate a relation of two masses Madm 
and My (Mk) to a moment equation of the relativistic Euler equation Eq. (A13). 

In the Newtonian theory, we usually check the accuracy of numerical solutions by the virial relation. Since the 
relation is not trivially satisfied in numerical solutions, violation of this relation can be used to estimate the magnitude 
of the numerical error of equilibria. Motivated by this idea, a virial relation is also derived for axisymmctric equilibrium 
states in general relativity [34] , and it is subsequently used to check accuracy of numerical solutions for rotating neutron 
stars [35]. The virial relation has been also derived for binary neutron stars in quasiequilibrium in conformally flat 
spacetimes [21] and applied for monitoring accuracy of numerical solutions in [36]. The virial relation, e.g., Eq. (All), 
derived here will be used when checking the accuracy of nonaxisymmetric numerical solutions. 
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